% This function looks for hot pixels in an image from sCMOS chip) and finds
% hot pixels and replaces their values with the mean value of the pixels
% around them

function [r,sortedHotpixelInd,numremoved] = removeHotPixels(im,thresh)
if nargin < 2 || isempty(thresh)
    thresh = .2;    
end

%remove some nonsense values:
bw = isnan(im);
im(bw) = median(im(:,1),'omitnan');

stdim = double(stdfilt(im));
medim = double(medfilt2(im));
stdim = stdim./medim;
im(abs(im)==inf) = medim(abs(im)==inf);
im(bw) = medim(bw);

temp = abs(stdim) > thresh;
SE = strel('diamond',1);
temp = imerode(temp,SE);
temp(abs(im)==inf) = 1;
numremoved = sum(temp(:));
hotpixelStdValues = stdim(temp);
sortedHotpixelInd = find(temp);
[~,I] = sort(hotpixelStdValues,'descend');
sortedHotpixelInd = sortedHotpixelInd(I);
r = im;
r(temp) = medim(temp);
